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: Abstract 

. The surface of a molecule determines much of its chemical and physical prop- 



erty, and is of great interest and importance. In this Letter, we introduce the 
concept of molecular multiresolution surfaces as a new paradigm of multiscale bi- 
ological modeling. Molecular multiresolution surfaces contain not only a family 
of molecular surfaces, corresponding to different probe radii, but also the solvent 
accessible surface and van der Waals surface as limiting cases. All the proposed 
■ surfaces are generated by a novel approach, the diffusion map of continuum sol- 

^ , vent over the van der Waals surface of a molecule. A new local spectral evolution 

^ kernel is introduced for the numerical integration of the diffusion equation in a 

single time step. 

Molecular surface [23] is one of the most important concepts in modern biochemistry 
and molecular biology. Albeit molecular surfaces are merely a conventional imitation 
of molecular boundaries, it has been shown that such models can often explain fun- 
damental physical and chemical effects. For example, the stability and solubility of 
macromolecules, such as protein, DNA and RNA, are determined by how their surfaces 
interact with solvent and/or other surrounding molecules. Therefore, the structure and 
function of macromolecules depend on the features of their molecular surfaces. Impor- 
tant applications of molecular surfaces [23], as well as solvent accessible surfaces [17] and 
van der Waals surfaces, are flaring in protein folding [20, 28], protein-protein interfaces 
[7], oral drug classification [2], DNA binding and bending [9], parameterization of heat 
capacity changes [20], macromolecular docking [16], enzyme catalysis [19], calculation of 
solvation energies [22], molecular dynamics [8], and implicit solvent models [1, 11, 26]. 

* Corresponding author. 
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Figure 1: Illustration of various surfaces. 

Molecular surface [23], also called solvent excluded surface, is the surface that en- 
closes the solvent excluded volume [5], and can be defined as a smooth envelope traced 
out by the surface of a probe sphere rolled over the molecular van der Waals surface. 
The latter is a collection of all the unburied atomic surfaces of the molecule, see Fig. 
1. The estimation and analysis of solvent accessible surfaces and molecular surfaces 
have attracted much attention in the past three decades after the first calculation of 
molecular surfaces by Greer and Bush [15]. A variety of methods have been proposed, 
including the Gauss-Bonnet theorem [5, 24], closed-form analytical method [13], space 
transformation [12], alpha shape theory [18], triangulation [6, 34], Cartesian grid based 
methods [26], contour-buildup algorithm [29], binary tree [27], and parallel methods [30]. 
Despit of much success, in general, calculation of molecular surfaces is both topologically 
and computationally challenging due to the multiple overlap of atomic spheres, overall 
irregular structure of a macromolecule, and self-intersecting singularities [14, 27]. The 
situation where the probe sphere is simultaneously tangent to four atoms poses another 
challenge [10]. Molecular surface calculations are the bottleneck in the molecular dy- 
namics of macromolecules, where the molecular surfaces are computed billions of times 
in the course of the simulation. 

About four decades ago, Benoi Mandelbrot posed the profound question of how one 
measures the length of a coastline in his paper entitled "How long is the coastline of 
Great Britain- Statistical self-similarity and fractional dimension" [21]. The answer 
depends on what scale the measurement is being made — one arrives at a different 
length by using a different scale. It is easy to envision that the molecular surface of a 
molecule depends the scale of the probe. As biological phenomena occur over a wide 
range of length scales, molecular mult iresolut ion surfaces ought to provide the scientific 
community with a multiscale framework to reveal complex biological processes across 
scales. Unfortunately, such a multiscale framework does not yet exist. 

The objective of the present work is to introduce the concept of molecular mul- 
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Figure 2: Multiscale surfaces marked up by flow contours, (a) The diffusion map of a 
fractal (Courtesy: The Mandelbrot Explorer Gallery); (b) The multiresolution surfaces 
of a water molecule. 

tiresolution surfaces both as a multiscale framework for biological modeling and as a 
unified description of the van der Waals surface, solvent accessible surface and molec- 
ular surface. The essential idea is similar to the multiresolution profiles of a fractal, 
generated by restricted diffusion, see Fig. 2. Use is made to continuum solvent and 
its curvature-controlled diffusion. The mathematical description of the solvent diffu- 
sion process is formulated by appropriate reaction-diffusion equations. All geometric 
features of a molecular surface, such as convex spherical patches, saddle-shaped pieces 
of tori, and concave reentrant faces, are naturally accounted by the geometric property 
of partial differential equations with appropriate initial-boundary conditions. A family 
of molecular multiresolution surfaces are found to correspond to different probe radii, 
and are generated in the continuous distribution of solvent density. The van der Waals 
surface is a special case in which the probe radius vanishes. A fast computational al- 
gorithm is introduced to generate molecular multiresolution surfaces in a single step of 
time evolution. 

Without the loss of generality, we assume that the fluid flow of continuum solvent has 
a density distribution p{r, t) at a space-time location (r, t), where r e and t e [0, oo). 
For an arbitrarily small volume Q, the motion of the solvent density is characterized by 
density flux J(r, t) through the surface of the volume, dfl. In this work, we assume that 
the flux is given by generalized Fick's law 

J(r, t) = -J2 D,{K{p), r, t) VV^r, t), (1) 
where Dq{K{p),r,t) are diffusion coefficients depending in the Gaussian curvature K 
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of solvent density isosurfaces, and high order terms (i.e., q > 1) describe the influence 
of inhomogencity in the solvent density and flux-flux correlation in the density flux. 
In particular, the term VV^ can be regarded as an energy flux operator. A similar 
fourth order flux term was employed in the Cahn-Hilliard equation [3] to describe the 
evolution of a conserved concentration held during phase separation. The total solvent 
in the volume is given by J^p{r,t)dr. The rate of change of the solvent in the vol- 
ume, f^p{r,t)dr, is balanced by the solvent out flux, — jQ^J-n{r)dS, and solvent 
production, J^P{K{p), p,r,t)dr 

[ ^p{r,t)dr = - [ J-niv)dS+ [ PiKip),p,r,t)dr (2) 
Jn Jon Jn 

where n(r) is the unit outer normal direction at r. By using the divergence theorem, we 
convert the surface integral into a volume integral and obtain 

= J2 ^■D,{K{p),r, VV^V(r, t) + P{K{p),p, r, t). (3) 

Q 

The diffusion coefficients Dq{K{p),r,t) are position dependent because of the possible 
presence of inhomogencity, such multivalent ions [25] , and their curvature-dependence is 
designed for the geometric control of solvent density isosurfaces in the final diffusion map. 
For example, the solvent density distribution near the molecular surface can be made to 

preserve the curvature of the probe molecule by taking Dq{K{p), r, t) ~ exp ^— ^ , 

where Kp is the Gaussian curvature of the probe and Kg are constants. The solvent 
production can be either negative or positive, and is designed for further geometric 
control of the molecular surface. It can also be used to add or remove atoms at specific 
locations. High order gradient-controlled diffusion equations were proposed by Wei 
[31] for efficient image processing. Anisotropic reaction-diffusion equations have been 
introduced for shock-capturing in the context of computational fluid dynamics [32] . 

In this work, we deflne the problem of flnding molecular multiresolution surfaces as 
an initial-boundary value problem governed by the new curvature-controlled diffusion 
equation (3) in a domain Q that is large enough to inchide all the surfaces of interest 
of a molecule. The Direchlet boundary condition is employed at domain boundary dQ, 
i.e., p{r,t) = po, re dD,. Initially, there is no presence of solvent in the domain, 
i.e., p(r, 0) =0, r E fl. We set Dq{K{p),T,t) = indiscriminately in all the atomic 
volumes Vi = |7rrf , where is the i^^ atom of the molecule. We choose the evolution time 
sufficiently long so that the solvent density can reach the center of the computational 
domain. The diffusion of the solvent toward to the domain center leads to a family of 
profiles of solvent density distribution, which in turn gives rise to a family of molecular 
multiresolution surfaces, corresponding to different solvent probe radii and different 
resolutions. Note that the van der Waals surface, dQvws, corresponds to Vp = 0, and 
is obtained at the level set p = 0. Here we define the midway molecular surface as 
the multiresolution molecular surface whose convex surfaces locate at the middle of 
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the van der Waals surface and the solvent accessible surface. Such a surface should 
be useful as it unclear whether the solvent accessible surface or the solvent excluded 
surface better describes hydration effects [16]. In general, the nonlinear terms in Eq. 

(3) can be utilized to provide large gradient near the desirable molecular surface. This 
approach directly results in a Cartesian grid representation of the molecular surface. 
Such a representation is important for many applications, particularly for implicit solvent 
modeling of electrostatics in macromolecules [26]. 

In general, it might not be convenient to solve the nonlinear reaction-diffusion equa- 
tion (3). Thus simplifications are desirable. For homogeneous solvent, it is sufficient to 
set Di a piecewise constant Di — D{r), and Dq = Q ii q 1. For simplicity, we also drop 
the production term in the rest of discussion. Therefore, we arrive at a much simpler 
problem of diffusion 

^^ = V-D(r)Vp(r,t), reO (4) 

p{r,t) = po, re dQ 
p(r,0) = 0, ren 

D{r) = 1, re n/Vvwv, 

D{r) = 0, re Vvwv, 

where Vvwv is the molecular van der Waals volume given by Vvwv — {JiVi- However, 
the resulting diffusion equation does not yet admit a closed-form solution due to the 
irregular atomic surface information in D{r). Fortunately, the numerical solution of Eq. 

(4) is quite robust. For example, both implicit and explicit approaches, as well as various 
spatial discretizations can be used to generate the solvent density profiles near the van 
der Waals surface. 

The van der Waals surface generated by this approach is exact in the sense that 
its resolution is only limited by the finite mesh size. However, the molecular surface, 
or the solvent excluded surface computed in this manner could be distorted because 
the different diffusion path length for each different part of a molecular. To eliminate 
this error of molecular anisotropy, we set the computational boundary to be the solvent 
accessible surface {dQ = dQsAs)- To this end, we define solvent- accessible radius of i^^ 
atom as the sum of its atomic radius and the probe radius, rf^^ = Vi + r-p. The volume 
associated with solvent-accessible radius of the i^^ atom is given by Vf"^^ = [rf"^^]^. 
We initialize the solvent density on the desirable solvent accessible surface by setting 
p(r, 0) = 0, Vr e IJ • Vf^^ and p(r, 0) = po in rest of the domain. A large probe radius 
Tp can be chosen to ensure a variety of molecular multiresolution surfaces. Note that 
solvent accessible surface marked out by the solvent density p(r, 0) is also exact. 

To speed up the solution of Eq. (4), we introduce a fast local spectral evolution kernel 
(LSEK) method [33], which is able to analytically solve a class of reaction-diffusion- 
convection equations = D(t)^ + C{t)-^ + P(t) p in a single time step. The 
LSEK is exact in time and is of controllable accuracy in space. Its complexity at a fixed 
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(c) (d) 

Figure 3: Molecular multiresolution surfaces of cyclohexane and their comparison with 
those obtained by the MSMS (Solid lines: present; Circles: MSMS). (a) The solvent 
excluded surface; (b) The cross section at a; = 0.6; (c) The cross section at y = 0.6; (d) 
The cross section at z = 0.6. 
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Figure 4: Molecular multiresolution surfaces of the cell division protein. Top: The 
van der Waals surfaces; Middle: The solvent excluded surface; Bottom: The midway 
molecular surface. 
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level of accuracy scales as 0{N), where N is the number of grid points. The solution of 
Eq. (4) at an arbitrary point of space-time, {xi, i/j, Zk,t), is obtained in a finite difference 
manner 

My 

P{Xi,yj,Zk,t) = XI 5Z 5Z ^h.,'TAlxhx,t) 

lx=-M^ ly=-My h=-M^ 

Khy,ay{lyhy,t)Kh,,aXhhz,t)p{xi - Ixhx.Vj - lyhy,Zk - lzhy,0), (5) 

where 2Ma + 1 are the stencil width for a — x,y, z, ha the grid spacings, and Kh^^{x, t) 
the LSEK [33] 




Here h2n{x) is the Hermite function defined by the Hermite polynomial -f^2n, i-e., h2n{x) = 
ex]){—x'^)H2n{x), Mh is the highest degree of Hermite polynomial used, a is window pa- 
rameter and CTj = + 2 D{s)ds. We set Mh = 88, Ma = 32 and Ua = 3.05/ia in this 
work. 

As a proof of principle, we present two examples. In the first example, cyclohexane 
(C6H12), a test case recommended by Sanner et al. [27], is considered to validate the 
proposed method for calculating the solvent excluded surface. The van der Waals radii 
of the carbon and hydrogen atoms are taken as 1.7 and 1.2 A, respectively. We choose a 
grid of 200 points in each dimension and set po = 100, r^ — 1.5 A. The solvent excluded 
surface is extracted at p(r, 12) = 0.04, see Fig. 3a. We also compute the solvent excluded 
surface of cyclohexane by using the MSMS algorithm [27]. A comparison of three cross 
sections generated from the present approach and those from the MSMS is depicted in 
Figs. 3b-d. There is an excellent agreement among these results. 

We next consider a cell division and cell wall biosynthesis protein [4]. It has 9772 
atoms in eight polymer chains and 1328 residues. We choose the same set of computa- 
tional parameters as those described in the preceding paragraph. Three multiresolution 
molecular surfaces, the van der Waals surface, the the solvent excluded surface, and the 
midway molecular surface of the cell division protein are depicted in Fig. 4. It is seen 
that the midway molecular surface captures the symmetry and other main features of 
the protein structure while avoids unnecessary details. 

In summary, we have introduced a new concept, molecular multiresolution surfaces, 
for the multiscale modeling of biomolecules. The proposed concept provides a unified 
description of the van der Waals surface, solvent accessible surface and solvent excluded 
surface. A novel approach based on the controlled diffusion of continuum solvent density 
is proposed to generate the multiresolution surfaces. A fast local spectral evolution ker- 
nel is introduced to integrate the diffusion equation in a single time step. The proposed 
method does not need to discriminate buried atomic surfaces from unburied ones during 
the computation. It by passes the difficulty of singularities in other existing methods. 
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Numerical experiments are carried out to validate and demonstrate the proposed ap- 
proach. The formulation and examples given here hopefully will spur the interest of 
other investigators in this and related fields. 
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